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Abstract 


A computationally efficient 3DHZETRN code with enhanced neutron and light ion (Z < 2) 
propagation was recently developed for complex, inhomogeneous shield geometry described by 
combinatorial objects. Comparisons were made between 3DHZETRN results and Monte Carlo 
(MC) simulations at locations within the combinatorial geometry, and it was shown that 
3DHZETRN agrees with the MC codes to the extent they agree with each other. In the present 
report, the 3DHZETRN code is extended to enable analysis in ray-trace geometry. This latest 
extension enables the code to be used within current engineering design practices utilizing fully 
detailed vehicle and habitat geometries. Through convergence testing, it is shown that fidelity in 
an actual shield geometry can be maintained in the discrete ray-trace description by 
systematically increasing the number of discrete rays used. It is also shown that this fidelity is 
carried into transport procedures and resulting exposure quantities without sacrificing 
computational efficiency. 

Introduction 

Primarily due to its high computational efficiency, the HZETRN code has been widely used throughout 
vehicle and habitat shield optimization activities, from preliminary architecture and trade studies through final 
design [Wilson et al. 2004a, b; Qualls et al. 2001; Shavers et al. 2004] and as one element of multidisciplinary 
optimization software [Wilson et al. 2003]. Further applications have been found in space environmental monitoring 
and validation [Hugger et al. 2003; Slaba et al. 2011a, 2013] and as a fundamental part of astronaut risk estimation 
for radiation safety monitoring [Cucinotta et al. 2008, 2013]. Despite the broad applicability, one limitation of the 
HZETRN model is the use of the straight-ahead approximation wherein all particles are assumed to travel along a 
common axis. This approximation is known to be inaccurate for low energy neutrons and light ions (Z < 2) produced 
at broad angles relative to the incoming primary radiation. 

Three-dimensional (3D) corrections to the HZETRN formalism for neutrons and light ions have recently 
been implemented, resulting in the 3DHZETRN code [Wilson et al. 2014, 201 5a, b]. These improved transport 
procedures have been compared to state of the art Monte Carlo (MC) codes, Geant4 [Agostinelli et al. 2003], 
FLUKA [Fasso et al. 2005, Battistoni et al. 2007] and PHITS [Sato et al. 2006, 2013], using both simplified 
spherical geometry [Wilson et al. 2014, 2015a] and more complex combinatorial geometry configurations [Wilson 
et al. 2015b]. In these benchmarks, it was generally concluded that 3DHZETRN agrees with the MC codes to the 
extent they agree with each other. 

The 3DHZETRN code currently allows the improved transport procedures to be evaluated within general 
combinatorial geometry that may be defined in terms of solid objects (spheres, boxes, cylinders, ellipsoids) with 
unlimited nesting (i.e. objects within other objects). This capability was tested against MC simulations by 
considering a cylindrical outer shield with six internal objects, including a tissue sphere [Wilson et al. 2015b]. 
Comparisons again showed that 3DHZETRN was in good agreement with the MC codes without sacrificing 
computational efficiency even in the more complex geometry considered. 

In principle, the current capabilities of 3DHZETRN would support the broad range of design studies and 
risk assessment applications discussed at the beginning of this section. However, using only a limited set of 
combinatorial solids to describe realistic and highly complex vehicle designs, such as the International Space 
Station, with thousands of internal parts would likely lead to unnecessary approximations and related errors in the 
geometric definitions. An alternative is to directly consider complex vehicle geometric descriptions typically 
developed in computer aided design (CAD) or finite element method (FEM) frameworks using ray-tracing methods. 
In such calculations, a ray-trace is performed over a large number of rays emanating from a target point within the 
geometry and extending outwards towards the exterior boundary. Along each ray, the length and type of each 
material traversed is stored. The resulting distribution of thicknesses is often referred to as a vehicle ray-trace or 
thickness distribution. Within the straight-ahead approximation of HZETRN, dosimetric quantities are quickly 
evaluated at the target point using either interpolation [Wilson et al. 1997] or ray-by-ray methods [Slaba et al. 
2011b] coupled to the ray-trace geometry. A major advantage of this approach is that shield design is allowed to 
evolve in a more suitable CAD or FEM environment, instead of in the generally more restricted geometric 
definitions allowed by either MC transport codes or 3DHZETRN. In this work, the 3DHZETRN code is extended to 
allow evaluation within geometries defined simply by a vehicle thickness distribution file, referred to herein as a 
ray-trace geometry definition. This extension allows 3DHZETRN to be easily interfaced with methods based on 
CAD or FEM models that are already in common engineering practice. No modifications to the transport formalism 
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are required for this extension of code capability. Only the ray-tracer used within 3DHZETRN along specific 
transport stream directions [Wilson et al. 2014] was modified. 

To test this new capability, the combinatorial geometry studied previously by Wilson et al. [2015b] will be 
replaced by the corresponding ray- trace geometry definition evaluated at different target point locations. 
Convergence tests are performed on the number of rays used in the thickness distribution to evaluate the associated 
geometric discretization error. It is shown that more precise geometric representations and transport code solutions 
are achieved as the number of rays in the thickness distribution is increased, and the current practice of using ~10 3 
rays may be sufficient in some cases. Results generated with 3DHZETRN from the ray-trace geometry are compared 
to 3DHZETRN results and MC simulations generated in the original combinatorial geometry [Wilson et al. 2015b]. 
It is found that neither transport code accuracy nor efficiency is sacrificed with the extension to ray-trace geometry. 
This latest transport code advancement directly enables the use of more complex geometric models so that simple 
mapping of the present methodology into more realistic applications is achieved. 

3D Transport through Ray-trace Geometry 

In this section, the ray-tracer used in 3DHZETRN to estimate material types and thicknesses traversed 
through a ray-trace geometry definition is described. The transport formalism previously described by Wilson et al. 
[2014, 2015a,b] is not altered. Fig. 1 shows a case in which the external radiation environment is assumed to 
impinge uniformly along the direction Qo onto the geometry, and the solution is being evaluated at the point x along 
the stream direction Q. Transport procedures require rapid evaluation of the materials and thicknesses traversed 
along an arbitrary direction extending from any point within the geometry to the outer boundary (along -Qo from xo 
in this example). Combinatorial geometry allows this evaluation to be performed analytically in closed form as 
described by Wilson et al. [2015b]. In ray- trace geometry, the exact geometry specifications are replaced with a 

discrete set of unit ray directions, r. , emanating from x , along which only the thickness, , and material types 

traversed by the ray are known. The superscript (j ) is used here to denote that in complex geometry, there may be 
many distinct material types, appearing along a given ray, i. 


Qo 



Fig. 1. General ray-trace geometry in which 3D transport procedures are evaluated. 


In order to efficiently evaluate traversed materials and thicknesses through ray-trace geometry, some 
requirements need to be defined. First, to enable rapid searching through the distribution of rays used to define the 
ray-trace geometry, it is required that these rays be generated according to a latitude-longitude grid, as shown in Fig. 
2. By using typical spherical coordinate transformations, the polar (longitude) and azimuthal (latitude) angles for 
each ray may be uniquely determined. The right side of Fig. 2 shows that these angles form a rectangular grid of 
points in two-dimensional space fully compatible with simple and efficient search algorithms used in bivariate 
interpolation. For an arbitrary ray-direction, its four nearest neighbors in the underlying ray distributions are easily 
found by simple search algorithms in two-dimensional space instead of computationally expensive comparisons 
against all ray directions in 3D space. 


2 



Fig. 2. Latitude-longitude distribution of rays mapped onto unit sphere (left) and onto two dimensional grid of polar 

and azimuthal angles (right). 


The only other requirement placed on the ray-trace evaluation mode is that at least one direction vector in 
the ray distribution must line up exactly with the incident direction of the incoming radiation environment, Qo. This 
requirement is implicitly required in existing engineering design practice already coupled to HZETRN in the 
straight-ahead approximation. The main reason for including this requirement in the ray-trace mode of 3DHZETRN 
is to ensure that evaluation of the forward flux component along Qo, which can dominate the overall exposure in 
many cases, is performed with precise geometric definition obtained directly from the native CAD or FEM 
geometry. 

After satisfying these computational requirements, the ray-trace procedure carried out within 3DHZETRN 
over the ray-trace geometry is relatively simple. Consider the problem of ray-tracing from a point x\ along a 
direction -Qo, as shown in Fig. 3. To begin, the ray direction, v\ (green dashed line in figure) connecting x\ to x is 
constructed and the polar and azimuthal angles for v\ are computed. The four nearest neighbors of v\ in the 
underlying ray distribution are then determined by quickly searching the 2D grid of values shown in Fig. 2. In the 
2D plot of Fig. 3, the visible neighbors are shown as blue solid lines. The closest of the neighbors, r c , is then 
selected, and the material type at a distance of d\ = |vi| along r c is obtained directly from the ray-trace geometry file. 
This material type is assigned to the point x\. A step of 0.1 mm is then taken in the direction of -Qo to define a new 
point, X 2 , and the process is repeated to determine the material type for point X 2 . 


J~2o 



This process is repeated until the distance d n = \v n \, where v n is the ray connecting x n to x, and x n is the point 
associated with step n in this procedure, is larger than the known distance along the nearest neighbor. This indicates 
that the outer boundary of the geometry has been reached. The end result of this procedure is a list of n material 
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types corresponding to points spaced 0.1 mm part in the geometry along the ray direction -Qo. Adjacent steps with 
identical material types are collapsed into a single step, so that the final information returned to the transport 
procedures is a list of contiguous material layers extending from the point of interest to the outer boundary in the 
specified direction. In the following section, various convergence tests are performed to evaluate functionality and 
accuracy of the described procedure. 


Inhomogeneous Geometry 


In this section, various convergence tests and comparisons are performed to test the ray-trace geometry 
mode of 3DHZETRN. In these tests, the combinatorial geometry of Wilson et al. [2015b], shown in Fig. 4, was ray- 
traced to generate the thickness distributions needed as input into the transport code. A precise definition of the 
original combinatorial geometry definition follows, but it should be understood that the analysis presented in this 
section is for transport code results using the ray-trace geometry definition derived from the combinatorial 
geometry. 

Fig. 4 shows a 30 cm diameter tissue equivalent sphere [ICRU 1993] placed within a 2.47 cm (5 g/cm 2 ) 
thick cylindrical aluminum shell (outer diameter 75 cm, outer length 175 cm) with 4.9 cm (10 g/cm 2 ) thick end caps. 
In addition to the International Commission on Radiation Units and Measurements (ICRU) sphere at the center of 
the aluminum cylindrical shell, an interior aluminum box with dimensions 46.3x14.8x14.8 cm (upper right corner in 
Fig. 4) and a second aluminum box with dimensions 22.2x22.2x37.0 cm (lower left comer in Fig. 4) are also placed 
inside the cylinder. The dimensions of these boxes in units of g/cm 2 are 125x40x40 and 60x60x100, respectively. 
Target points (detectors) are placed within the ICRU tissue sphere, similar to the studies of Wilson et al. [2014, 
2015a]. The detector locations appearing at the top and bottom of the ICRU tissue sphere will be referred to in 
subsequent figures, tables, and discussion as z = 0 and z = 30, respectively (z is parallel to Qo). Additional target 
points within the ICRU tissue sphere are therefore defined by their location between the top and bottom detectors. 
For example, the detector location in the center of the sphere is referred to as z = 15. 

For all calculations, the external radiation environment (external source boundary condition) is positioned 
above the geometry and directed uniformly downward along Qo as indicated by the arrows in Fig. 4. The boundary 
condition, or external source, is the Webber representation of the February 23, 1956 solar particle event (SPE) 
[Webber 1966]. The Webber 1956 SPE is given by an exponential momentum spectmm with momentum parameter 
po = 100 MV and 10 9 protons/cm 2 above 30 MeV. 

To begin the process of testing transport procedures in ray-trace geometry, the number of rays required to 
adequately represent the geometry is examined. It is important to note that the convergence tests and related 
discussion and conclusions in this work are specific to the geometry being considered. More complex geometry 
would necessitate similar testing procedures to quantify the discretization error associated with representing the 
geometry with a finite number of rays. 
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A - cylindrical aluminum shell 
B - ICRU tissue sphere 
C - aluminum box 
D - aluminum box 
E - detectors (target points) 




Fig. 4. Complex geometry used in the current study: ICRU tissue sphere shielded by an aluminum cylindrical shell 

with internal parts. 
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Fig. 5 gives a sequence of images showing converging geometric resolution with increasing number of 
rays, N my , in the original ray-trace file. These images were created by first generating a series of ray-trace files, each 
one corresponding to a target point placed at the center of the ICRU tissue sphere and with a different N ray value. 
These files were then individually input into 3DHZETRN, and a large number of target points and ray directions 
were arbitrarily chosen. Each target point and direction was passed to the ray-tracer described in the previous section 
to evaluate approximate location of material interfaces. The location of the material interfaces, as determined from 
the internal ray- tracer, were then plotted as scatter points to create the images in Fig. 5. 



Fig. 5. Geometric representation associated with ray-trace geometry (N my = 970, 3874, 15490, 61954) and 

combinatorial geometry. 

Two observations should be made in connection to Fig. 5. First, objects closer to the original target point at 
the center of the ICRU sphere at which ray-traces were evaluated are in clear focus even with a small number of 
rays. As the distance from the original target point increases, so too does the lateral distance between adjacent rays, 
leading to greater geometric distortion. If the target point used in the original ray-trace files were placed in the upper 
right box, for example, greater distortion would occur in the ICRU sphere, and the upper right box would be in 
clearer focus regardless of the chosen N my value. Second, although considerable geometric distortion occurs for the 
lower number of rays, the geometric representation is seen to systematically focus at the larger number of rays, 
indicating that a converging sequence of transport solutions should be expected for increasing N my values. 

An important point to keep in mind is that the geometric distortion shown in Fig. 5 only impacts the 
isotropic components of the neutron and light ion spectra. By requirement, as noted in the previous section, the 
forward component of the transport solution is evaluated directly over a ray for which exact material and thickness 
information from the original geometry is known, and no approximate ray-trace methods are utilized. In many cases, 
the forward component of the flux dominates the overall exposure, so that total exposure quantities and fluence 
spectra would be expected to be reasonably accurate. Moreover, geometry distortion far from the target point of 
interest, as shown in Fig. 5, will necessarily lead to errors in the spatial distribution of isotropic neutron source terms 
driving 3D transport procedures. This can be seen in Fig. 6, where the integral isotropic neutron source induced 
within the geometry by the Webber 1956 SPE is plotted. The source terms far from the evaluation point play less of 
a role in determining the isotropic neutron and light fluence spectra than the source terms very near the evaluation 
point, for which geometric distortion is minimal. This also indicates that the distortion shown in Figs. 5 and 6 may 
not cause substantial error in the induced particle spectra within the geometry. Still, the impact of this apparent 
geometric distortion on the isotropic components of the neutron and light ion spectra needs to be further 
investigated. 


5 


A U =970 r-* 



Fig. 6. Integral ( E > 1 MeV) isotropic neutron source (particles/(g-event)) induced by the Webber SPE using ray- 
trace geometry ( N my = 970, 3874, 15490, 61954) and combinatorial geometry. 


The effect of geometry discretization on radiation field quantities is now examined. Wilson et al. [2015b] 
previously established that N= 26 transport streams are required to evaluate the radiation fields in the combinatorial 
geometry being considered presently. Although results are not shown here, it has been verified that N = 26 transport 
streams are adequate in the corresponding ray-trace geometry as well, and results given in subsequent tables and 
figures were all generated with N = 26 transport stream directions. Of greater interest is the number of rays required 
in the ray-trace geometry file to adequately describe the actual geometry. Figs. 5 and 6 have already provided some 
indication of preliminary requirements, and now dosimetric quantities and fluence spectra are evaluated to provide 
quantitative assessments. 

Fluence spectra evaluated for various N ray values are shown in Fig. 7, and dose and dose equivalent results 
are given in Tables 1 and 2. Despite the geometry and neutron source distortions seen in Figs. 5 and 6, fluence 
spectra are nearly identical with results appearing as overlapping on the log-log scale plot of Fig. 7. 




Fig. 7. Fluence spectra evaluated at the top (left) and bottom (right) of the ICRU sphere for both combinatorial 

representation and ray-trace geometry approximation. 
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Table 1. Dose at points in tissue sphere for various N ray values used in ray-trace geometry. Results from Wilson et al. 
[2015b] for the combinatorial geometry are also provided for comparison. 


z (g/cm 2 ) 


Dose (cGy/event) 


N mv = 970 

Nray = 3874 

Nray = 15490 

Nray = 61954 

Combinatorial 

0 

42.78 

42.78 

42.79 

42.79 

42.79 

0.007 

42.48 

42.48 

42.48 

42.48 

42.49 

0.3 

38.74 

38.74 

38.75 

38.75 

38.75 

1 

32.38 

32.38 

32.39 

32.39 

32.39 

5 

14.59 

14.60 

14.60 

14.60 

14.59 

25 

1.55 

1.55 

1.55 

1.55 

1.55 

29 

1.13 

1.13 

1.13 

1.13 

1.13 

29.7 

1.08 

1.08 

1.08 

1.08 

1.07 

29.993 

1.05 

1.05 

1.05 

1.05 

1.05 

30 

1.04 

1.04 

1.04 

1.04 

1.04 


Table 2. Dose equivalent at points in tissue sphere for various N my values used in ray- trace geometry. Results from 
Wilson et al. [2015b] for the combinatorial geometry are also provided for comparison. 


Dose Equivalent (cSv/event) 


z ) 

Nray = 970 

Nray = 3874 

N mv = 15490 

Nray = 61954 

Combinatorial 

0 

65.89 

65.90 

65.94 

65.97 

65.99 

0.007 

62.39 

62.40 

62.45 

62.48 

62.51 

0.3 

55.22 

55.24 

55.30 

55.33 

55.36 

1 

46.03 

46.04 

46.10 

46.12 

46.15 

5 

20.78 

20.81 

20.82 

20.82 

20.81 

25 

2.62 

2.61 

2.60 

2.59 

2.59 

29 

2.08 

2.07 

2.06 

2.06 

2.04 

29.7 

2.03 

2.01 

2.00 

2.00 

1.98 

29.993 

1.99 

1.97 

1.96 

1.96 

1.94 

30 

1.90 

1.89 

1.87 

1.87 

1.85 


Slight discrepancies can be seen in the neutron spectra below 10 MeV/n and are quantified in Fig. 8 where 
the relative difference between the neutron spectra obtained in the combinatorial geometry and various ray-trace 
geometries are plotted. Relative differences for other light ions were found to be smaller than what is shown in Fig. 
8 for neutrons. It can be seen that relative differences are within +20% across the full energy range and diminish 
with increasing N ray values. Relative differences are also reduced at the highest energies where the forward 
component of the flux dominates the solution, and geometric distortion shown in Fig. 6 has negligible impact. 

In Table 1, it is found that dose converges quickly with increasing N my , and only negligible differences are 
found between results derived from ray- trace and combinatorial geometry descriptions. Dose equivalent values 
derived from the ray-trace geometry are only a few percent different than the values obtained within the more 
detailed combinatorial geometry. The slightly larger variation observed for dose equivalent is attributed to errors in 
the low energy light ion fluence spectra as shown in Fig. 7. 

These errors can be more clearly seen in Fig. 9, where the alpha fluence from Fig. 7 has been shown on an 
expanded scale. The low energy alpha fluence at z = 30 shows very slight systematic errors as a function of the N my 
value. It can be seen that the largest N ray value of 61954 yields results almost identical (within 2%) to the 
combinatorial geometry results, suggesting that even small-scale detail in the particle fluence spectra found at 
evaluation points near material interfaces can be captured within the ray-trace geometry, provided that enough rays 
are used. 
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Fig. 9. Alpha fluence from Fig. 7 with vertical axis in linear scale to improve clarity. 


It is encouraging that the differences associated with varying N ray show little effect on the spectra, dose, and 
dose equivalent and that computational efficiency has been maintained in the ray-trace geometry mode, as indicated 
in Table 3. Fig. 10 gives a final comparison between MC results and 3DHZETRN in combinatorial mode [Wilson et 
al. 2015b] along with new results in ray-trace mode with N ray = 61954. The ray-trace derived results are clearly in 
good agreement with the MC simulations and with the combinatorial results. It can also be seen that variation 
amongst the MC codes, associated almost entirely with differing nuclear interaction models, greatly exceeds the 
variation associated with various N my values used in 3DHZETRN in ray-trace mode, as shown in Fig. 7. This 
indicates that the geometric distortion associated with ray-trace geometry definitions leads to transport model errors 
that are much smaller than current uncertainties in the underlying nuclear physics models used in the codes. 

Table 3. Total CPU time required for transport calculations. 


Code 

Combinatorial 

Nray = 970 

Nray = 3874 Nray = 15490 

Nray = 6 1954 

CPU time (sec) 

59 

66 

68 71 

83 
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Fig. 10. Fluence spectra evaluated at the top (left) and bottom (right) of the ICRU sphere for both combinatorial 
representation and ray-trace geometry approximation compared to MC simulations. 


Summary 

In this work, extensions to the 3DHZETRN code were presented, allowing the code to be evaluated in ray- 
trace geometry. This latest update enables the code to be coupled to highly detailed vehicle and habitat designs that 
are commonly developed, maintained, and optimized in CAD or FEM environments. All that is needed for the 
coupling between 3DHZETRN and fully detailed geometry is the generation of a single thickness distribution 
through ray-tracing procedures already implemented in current engineering design practice. Minimal requirements 
are imposed to enable efficient coupling and were described. These requirements allow efficient look-up procedures 
to be carried out within 3DHZETRN over the ray- trace geometry and are actually consistent with current 
requirements already in use with HZETRN in the straight-ahead approximation. 

To test the updated transport code, the combinatorial geometry of a prior study [Wilson et al. 2015b] was 
considered, and convergence tests were performed. These tests revealed that increasing the number of rays used in 
the underlying ray-trace geometry description sufficiently reduced geometric discretization error and spatial 
distortion of the isotropic neutron source term. Transported exposure quantities such as dose, dose equivalent, and 
fluence were also examined to quantify error associated with reducing precise combinatorial geometry to discrete 
ray-trace representations. It was shown that both dose and dose equivalent values can be precisely calculated in the 
ray-trace geometry with negligible error compared to values derived from the exact combinatorial geometry. 
Comparisons to previous MC simulations [Wilson et al. 2015b] also showed that uncertainties in the nuclear physics 
models greatly exceed the error associated with discrete ray-trace geometric representations, indicating that the 
updated code is suitable for a broad range of engineering applications at least in the context of present nuclear 
interaction knowledge. Clearly, an increase in fidelity of nuclear models and their validation with laboratory 
experiments are needed. 
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